package scu.maqiang.homogenization;

import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.fes.FES2T61;
import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.Mesh2T6;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class PossionEigenValue2DT3 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
	    int N = 16;
	    double ep= 1.0 / N;
		
		Mesh2T3 cellMesh = new Mesh2T3().square(32, 32);
		ScalarFunc regionCellFunc = (xy, llable, param) -> {
            return (Math.abs(xy[0] - 0.5) < 0.25 && Math.abs(xy[1] - 0.5) < 0.25)? 1.0:-1.0;
        };
	    cellMesh.setDomainLabel(regionCellFunc, null, 10);
	    FES2T31 fsCell = new FES2T31(cellMesh);
	    double[][] FOCS = new double[2][fsCell.GetNdof()];
		double[][] SOCS = new double[4][fsCell.GetNdof()];
		double[][] D = new double[2][2];
		SOTS2DT361 ss = new SOTS2DT361(fsCell, null);
		double[][] coef1 = {{0.01, 0}, {0, 0.01}};
		double[][] coef2 = {{1.0, 0}, {0, 1.0}};
		MatrixFunc mf = (xyz, label, param) -> label == 10? coef1 : coef2;
		ScalarFunc rhoF = (xyz, label, param) -> label == 10? 0.5 : 1.0;
		double[] homoRho = {0.0};
		ss.EigenValueCellSolution(mf, new ScalarFunc[] {rhoF}, null, FOCS, D, homoRho, SOCS);
		
		System.out.println("min: " + MVO.min(FOCS[0]) + "\tmax: " + MVO.max(FOCS[0]));
		System.out.println("min: " + MVO.min(FOCS[1]) + "\tmax: " + MVO.max(FOCS[1]));
		System.out.println("Homogenized Kxx: ");
		System.out.println(MVO.toString(D));
		System.out.println("Homogenized Density");
		System.out.println(homoRho[0]);
		System.out.println("min: " + MVO.min(SOCS[0]) + "\tmax: " + MVO.max(SOCS[0]));
		System.out.println("min: " + MVO.min(SOCS[1]) + "\tmax: " + MVO.max(SOCS[1]));
		System.out.println("min: " + MVO.min(SOCS[2]) + "\tmax: " + MVO.max(SOCS[2]));
		System.out.println("min: " + MVO.min(SOCS[3]) + "\tmax: " + MVO.max(SOCS[3]));
		cellMesh.toTecplot("EigenValueCellFOCS.dat", FOCS);
		cellMesh.toTecplot("EigenValueCellSOCS.dat", SOCS);

		Mesh2T6 homoMesh = new Mesh2T6(cellMesh);
		FES2T61 homofs = new FES2T61(homoMesh);
		SRMatrix homoA = new SRMatrix(homofs.GetNdof());
		SRMatrix homoM = new SRMatrix(homofs.GetNdof());
		
		int eigenNumber = 10;
		double[] homoD = new double[eigenNumber];
		double[][] homoV = new double[eigenNumber][homofs.GetNdof()];
		double[] homoX = new double[homofs.GetNdof()];
		Arrays.fill(homoX, 1.0);
		MatrixFunc homoC = (xy, label, param) -> D;
//		homofs.assembleHeatStiff(new double[] {D[0][0]}, BVPType.COMMON, homoA);
		homofs.assembleStiff(homoC, null, BVPType.CONSTITUTE_MATRIX_COMMON, homoA);
		homofs.assembleMass(homoRho, BVPType.CONSIST_MASS, homoM);
		homofs.applyBCZero(homoA, homoX, Direct.X, BCType.RCE, 1, 2, 3, 4);
		EigenSSolver solver = new EigenSSolver();
		solver.inversePowerTran_PCGSSOR(homoA, homoM, homoX, 0.0, eigenNumber, homoD, homoV);
		System.out.println("Homo EigenValues: ");
		System.out.println(MVO.toString(homoD));
		homoMesh.toTecplot("HomoEigenVectors.dat", homoV);
		
		
		Mesh2T3 fineMesh = new Mesh2T3().square(512, 512);
		ScalarFunc regionFineFunc = (xy, rrlabel, param) -> {
            double xCell = xy[0] * N - (int)(xy[0] * N);
            double yCell = xy[1] * N - (int)(xy[1] * N);
            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
        };
        fineMesh.setDomainLabel(regionFineFunc, null, 10);
        
		FES2T31 finefs = new FES2T31(fineMesh);
		SRMatrix fineA = new SRMatrix(finefs.GetNdof());
		SRMatrix fineM = new SRMatrix(finefs.GetNdof());
//		
		double[] fineD = new double[eigenNumber];
		double[][] fineV = new double[eigenNumber][finefs.GetNdof()];
		double[] fineX = new double[finefs.GetNdof()];
		Arrays.fill(fineX, 1.0);
		ScalarFunc cf = (xyz, label, param) -> label == 10? 0.01: 1.0;
		finefs.assembleStiff(new ScalarFunc[] {cf}, null, BVPType.COMMON, fineA);
//		A.sort();
//		System.out.println(A);
		finefs.assembleMass(new ScalarFunc[] {rhoF}, null, BVPType.CONSIST_MASS, fineM);
		//M.sort();
		//System.out.println(M);
		finefs.applyBCZero(fineA, fineX, Direct.X, BCType.RCE, 1, 2, 3, 4);
		//System.out.println(A);
		solver.inversePowerTran_PCGSSOR(fineA, fineM, fineX, 0.0, eigenNumber, fineD, fineV);
		//solver.inversePowerTran_PLU(A, M, fs.getnBand(), x, 2.0, N, D, V);
		System.out.println(MVO.toString(fineD));
		fineMesh.toTecplot("FineEigenVectors.dat", fineV);
//		
		double[][] V0 = new double[eigenNumber][finefs.GetNdof()];
		double[][] V1 = new double[eigenNumber][finefs.GetNdof()];
		double[][] V2 = new double[eigenNumber][finefs.GetNdof()];
		for (int i = 0; i < eigenNumber; i++) {
			ss.assembleHeatSOTSSolution(finefs, V0[i], V1[i], V2[i], homofs, homoV[i], fsCell, FOCS, SOCS, N, ep);
		}
//		
		
		double[] D1 = new double[eigenNumber];
		double[] D2 = new double[eigenNumber];
		ss.assembleEigenValueSOTSSolution(finefs, V0, V1, V2, rhoF, homoD, D1, D2, ep);

		fineMesh.toTecplot("HomoEigenVectors.dat", V0);
		fineMesh.toTecplot("FotsEigenVectors.dat", V1);
		fineMesh.toTecplot("SotsEigenVectors.dat", V2);
		
		System.out.println("Fine          Homo          FOTS          SOTS");
		for (int i = 0; i < eigenNumber; i++) {
			System.out.println(fineD[i] + "\t" + homoD[i] + "\t" + D1[i] + "\t" + D2[i]);
		}
		
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");

	}

}
